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ABSTRACT 

We present results of modeling the broadband SED and multiwavelength variability of the 
brigh t FSRQ PKS 1510— 089 with our time-dependent multizone Monte Carlo/Fokker- Planck 
code dChen et al.ll201 lbl) . As the primary source of seed photons for inverse Compton scatter- 
ing, we consider radiation from the broad line region (BLR), from the hot dust of the molecular 
torus, and the local synchrotron radiation (synchrotron-self-Compton, SSC). We evaluate the 
viability of different Compton models by comparing simulated multiwavelength light curves 
and SEDs with one of the best observed flares by PKS 1510-089, in March 2009. The time- 
dependence of our code and its correct handling of light travel time effects allow us to fully 
take into account the effect of the finite size of the active region, and in turn to fully exploit the 
information carried by time resolved observed SEDs that are becoming increasingly available 
since the launch of Fermi. We confirm that the spectrum adopted for the external radiation 
field has an important impact on the modeling of the SED, in particular for the lower energy 
end of the Compton component which is observed in the X-ray band, which in turn is one of 
the most critical bands to assess the differences between EC and SSC emission. In the con- 
text of the scenario presented in this paper, where the flaring is caused by the increase of the 
number of relativistic electrons ascribed to the effect of the interaction of a portion of the jet 
(blob) with a shock, we can not firmly discriminate the three main scenarios for 7-ray emis- 
sion. However, results show clearly the differences produced by a more realistic treatment of 
the emitting source in the shape of SEDs and their time variability over relevant, observable 
time-scales, and demonstrate the crucial importance of time-dependent multi-zone models to 
advance our understanding of the physics of these sources, by taking full advantage of the 
wealth of information offered by the high quality data of current multiwavelength campaigns. 
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1 INTRODUCTION 

The spectral energy distributions (SED) of blazars usually show 
two major non-thermal components. The low energy one, peak- 
ing in the IR - optical - X-ray range is identified as synchrotron 
radiation. The origin of the high energy one, peaking in the 7- 
ray energy range (MeV to TeV) is less clear. Proposed ideas in- 
clude leptonic models, based on inverse Compton (IC) scatter- 
ing by the same electrons emitting the synchrotron radiation, and 
hadronic models, in which protons play a critical role in produc- 
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ing the high energy emission (Mannheim 1998; Rachen 2000; 
Sikora & Madeis kTl200ll; lArbeiter et al.l 120051 ; iLevinsonl l2006t 
Bottcher! |2007| ; iDermer & Lottl 1201 ll) . For the leptonic IC-based 
models, several sources of the target photons are possible and the 
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debate about which one is dominant and in which type of object 
has recently been reignited. If the seed photons are provided by 
the synchrotron radiation emitted at lower energy by the same IC- 
scattering electrons, it is referred to as synchrotron self -Compton 
(SSC). Scenarios in which the dominant contribution to the seed 
soft photon field for IC is provided by radiation emitted else- 
where are referred to as external Compton (EC) models. External 
sources of seed photons may include the photons from accretion 
disc, hot X-ray-emitting corona, broad emission line region, in- 
frared torus, host galaxy bul ge, and cosmic background radiation 
dGhisellini & Tavecchii 



xy bulge 
120091 



There are two major classes of blazars: BL Lac objects, which 
have featureless optical spectra, and Flat Spectrum Radio Quasars 
(FSRQ), which exhibit broad quasar-like emission lines. The pres- 
ence of these latter in FSRQs suggests that their jets are in an envi- 
ronment with a stronger external radiation field. Furthermore, de- 
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pending on the relative location of these sources external to the 
jet and of active jet region (blob), the emission from the exter- 
nal sources would be relativistically beamed and enhanced in the 
frame of the blob, possibly making them dominant over the locally 
produced synchrotron emission. Therefore, the EC model is fre- 
quently invoked to explain the emission of FSRQs dDermer et alj 
ll992tl'sikora et al .11 1994 1200^) . 

A major defining feature of blazars is their rapid and large 
variability, observed over the entire range of their continuum emis- 
sion (radio to TeV 7-rays). Simultaneous multiwavelength observa- 
tions and the correlation analysis of the observed multiwavelength 
variability can provide insights to the physics of particle processes 
and radiation mechanisms in the jet. The detailed observation and 
modeling of such variability has been performed ex tensively for 



High energy peaked BL L acs (HBL) such as Mrk 421 (Fossati et al 
20081; IChen et alj 1201 lbl) and PKS 2155-304 dAharonian et~ 
20071 ; Katarzvnski et al. 20081) , In the case of HBLs, such stud- 



ies are facilitated by their synchrotron SED peak falling in the X- 
ray range which can be observed with multiple X-ray satellites, 
while their high energy SED peak in TeV 7-ray is covered by 
ground based Atmospheri c Cherenkov Telescopes ( for a review, 
Hinton & Hofmann 2 0091; exa mples for Mrk 421, iFossati et alj 
2008l ; lAbramowski et al.ll2012l) . 

The launch of Fermi has re-opened the GeV 7-ray sky with un- 
precedented sensitivity and daily coverage. This energy band cov- 
ers a highly variable part of the SED right above the peak of the 
high energy compone nt of the SED of several bright FSRQ s, such 
as PK S 1510-089 dAbdo et aljfcoiOal) and 3C454.3 dAbdo et alj 
2009). Simultaneous coverage in other wavelength such as optical 
and X-rays provided us a chance to obtain multi-epoch SEDs and 
cross-band correlations. A deeper understanding of these time se- 
ries data sets requires time-dependent modeling with all light travel 
time effects (LTTEs) taken into account. 

The imp ortance of the LTTE in the stu dy of blazars has long 
been realized ( Chiab erge & Ghisellinil 1999b . The observed change 
of flux level on time-scale of hours in some blazars indicates 
that these effects must have a large impact on the variability of 
blazars. There have bee n some efforts to include these effects in 
the modeling of blaz ars (Kataoka et al. 2000; Sokolov et al. 2004; 



Sokolov & Marscher! 120051 ; iKatarzvfiski et al] 120081: iGraff et alj 
2008), but none of them have taken into account all of these effects 



in a generic 2D geometry. 

We have developed a time-dependent multizone code using 
the Monte Carlo method for radiation transport and the Fokker- 
Planck equation for electron evolution, which we first applied to a 
study of the correlated X-ra y/7-ray va riability of the HBL Mrk 421 
using a pure SSC model dChen et aljkoilbl) . In this paper we ex- 
tend the model to include external sources of IC seed photons and 
apply it to study the multiwavelength variability of an archetypical 
powerful FSRQ, PKS 1510-089. 



monitoring effort, complemented by more intense campaigns mo- 
tivated by flaring phases, has lead to the o bservation of sever al 
large correlated flares for PKS 1510-089 dAbdo et al]|2010bllal; 
ID Ammando et alj I20T1I ; iMarscher et al] 120101; iKataoka et al] 
l2008l ; |Pucella et alj|2008l) . 

As reference data for our simulation we choose the obser- 
vations of a high s t ate observed in 200 8-2009 and presented by 
lAbdo et al] d2010ah . lD ; Ammando et al] d201lh and lMarscher et al] 
d2O10l) . One particular flare at the end of March 2009 (peaking 
around March 25th, MJD 54917) is chosen as the benchmark for 
this study. 

We aim to reproduce several observational features by match- 
ing both the simulated light curves and SED with the observed 
ones. These features include: 



• Clear presence of BBB emission with the general two non- 
thermal continuum components SED. The BBB is evident in both 
the high and the low states of the jet non-thermal continuum emis- 
sion. 

• The time-scale of the flares, which were typically about 4 days 
(300 ks) at all wavelengths. 

• Infrared and 7-ray light curves show the stronger variations 
among the observed bands, with similar amplitude, up to a factor 
of 10. In the March 2009 flare, the infrared (R band) and 7-ray 
(FermifLAT) fluxes were strongly correlated, with no significant 
lags. 

• The variations in the SwiftfLWOT bands were less prominent 
than those in the infrared - optical bands. The optical/UV spectral 
shape became softer when the source brightness increased, consis- 
tent with the combination of a variable softer broad band continuum 
(synchrotron) and a non variable component peaked in the UV band 
(BBB). 

• The variability in the X-ray band luminosity was modest, 
within a factor of 2 over a period of 4 months en compassing 
the fla re that we selected for this study (Swiftf KRT in lAbdo et al] 
l2010al RossiXTEfPCA in lMarscher et al.l2010h . The spectrum was 
always very hard, with energy index ax < 0.6 (F v oc v~ a ) with 
very little variability. D uring the period aro und the benchmark flare 
it remained a x ^ 0.5 dAbdo et all2010al) . 

• The > 0.2 GeV spectral shape as measured by FermifLAT did 
not vary significantly through large lu minosity changes, r emaining 
around 07 ~ 1.5 (for a power law fit) (Abdo et aL2 010al) . 

These characteristics only coarsely summarize the true richness of 
information provided by the full multiwavelength and multi-epoch 
data set which thanks to our simulations we can try to exploit more 
deeply. Nevertheless, because they constitute a quicker and easier 
way of guiding the setup and evaluation of the simulations and we 
will refer to them when discussing the comparison of our simula- 
tion results with observations in the following sections. 



2 PKS 1510 089 

PKS 1510-089 is a FSRQ at a redshift of 2 = 0.361 
dThompsonetalJl99d) . It is one of the brightest and most variable 
sources detected by FermifLAT. A feature that can be interpreted 
as disk emission (big blue bump, BBB) is clearly visible in its op- 
tical/UV spectrum. VLBI observations of its jet show superluminal 
motio n with apparent speed up to 45c ((r) = 36. Ijorstad et al] 
l2005h . 

Since the advent of Fermi, the long-term multiwavelength 



3 SIMULATIONS 

3.1 Basic setup and model parameters 

Details and technical aspects of our Monte Carlo/Fo kker-Planck 
code are described in IChen et aT] d201 lbl IChenll2012i Chen et al. 
in preparation). The code uses the Monte Carlo method to track 
the production, travel, and Compton scattering of photons, while it 
solves the isotropic Fokker-Planck equation to follow the evolution 
of electrons. The major strength and unique feature of this code is 
that it takes into account all the LTTEs, internal to the source vol- 
ume (e.g. important for IC emission) and external, i.e. their effect 
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Table 1. Summary of model parameters 



Parameter 


n f / tol r 


blr 1 5 


bl r 1 5highgmm 


bl r 2 5 


torus 1 5 


SSC 


Jet bulk Lorentz factor, V 


1 J.U 


15.0 




25.0 


15.0 


10.0 


Zj \l\J CITIJ 


o.u 


o.u 


O.U 


1 J. JJ 


O.U 


J.U 


_rt cm ) 


o.u 


o.u 


o.u 


1U.U 


O.U 


1 7^ 

j. / j 


Magnetic field, B 


U. J 




U.J 


n 1 a 
u. lO 


U.Z 


U. 1 


Particle density (initial), n c (10 4 cm ^) 


z.oo 


Z.00 


Z.OO 


U. 14 


/ . J / 1U 


U.Ul 


L al L1L/1C CaUtlJJC L1111C ■aLdlC, CSC V / / 


1 


1 


1 


0.1 


015 


03 


PiTti r*l ( Hi "f~H i c f*A nr*r*f*l pnti r^n tim^-cr'nlf 1 ^ ^ / /^l 
L al L1L-1C ^UlllUSCy u.LLC1C1u.L1*J11 L1111C f»Ud.lC, ^^(^q \^ j *-■/ 


0.55 


55 


55 


55 


09 


19 


Electron pick up rate Qpjck (cm - '"' s ^~ ) 


0.1 


0.1 


0.1 


3.2 \0~ 6 


0.0191 


0.002 


Pick-up electrons energy, 7 p i c k 


5.0 


5.0 


5.0 


5.0 


50.0 


1200.0 


Shock injection: 7 m in,inj 




30 


90 


6 


300 


2000 


Shock injection: 7 max ,inj 




2 10 4 


2 10 4 


4 10 3 


2 10 5 


10 5 


Shock injection: power-law slope, pi n j 




3.2 


3.2 


3.2 


3.2 


3.2 


Shock injection: rate: Lj n - (10 44 erg s" 1 ) 




3.5 


2.0 


2.8 


5.0 


8.0 


^blr or R m (lQ l * cm) 


0.8 


0.8 


0.8 


0.8 


7.8 




/blr or /ir 


0.013 


0.013 


0.013 


0.0015 


0.5 





on the observed radiation which we will receive at different times 
depending on where it was emitted in the source. 

We model a jet active region (blob) as a cylindrical volume 
crossing a standing 'shock' as illustrated in FigQ] In the blob 
rest frame, where all calculations are performed, the shock moves 
through the cylindrical region with a speed equal to the bulk veloc- 
ity of the blob Vbuik ~ c. The cylindrical volume is divided evenly 
into zones in the radial and vertical directions (r and z coordinates, 
n r , n z ). In all runs presented in this paper, n r — 9 and n z = 30. 

At this stage the meaning of this shock is simply that of an 
agent affecting the properties of the simulation zone where it is at a 
given time. In the simulations presented in this paper it affects the 
electron distribution, namely it injects high energy particles into the 
zones it currently resides in. It does not affect the magnetic field, 
which is assumed to be homogeneous throughout the volume and 



Electron Density »t. 
Electron Spectrum ^ ' ' 

Magnetic Field B(r,z) 




shock, stationary 
in the BH frame 



bulk motion direction 
in the BH frame 



external radiation 



Figure 1. The geometry of the model. The volume is divided in the r and z 
directions in zones with their own electron distribution and magnetic field. 
We also schematically show the setup for the variability of the simulations 
with a shock. The hatched layer represents a stationary shock. The blob 
moves downward and crosses the shock front. For illustration purposes, we 
plot in lighter color shade the zones that crossed the shock at earlier times 
and have had some time to radiate the newly injected energy. In this repre- 
sentation the photons from the external radiation fields, beamed in the blob 
rest frame, enter the blob from the bottom surface. 



non varying. The shock has not thickness (however, our resolution 
is limited by dz) and is treated as a surface perpendicular to the z 

direction. 

W ith respect to the simulations presented in dChen et at] 
l2011bh , we added a few features to the code. The most relevant 
ones are the inclusion of a spatially diff use stochastic acce l eratio n 
process, similar to the one discussed by iKatarzvnski et ail ( feOOrJ) . 
of a particle escape ter m and of a particle 'pick-up' term (see also 



Tramacere et al. 



201 lh. The s e deve l opmen t s and their motivation 



Chen et al.l f201 lal) . IChenl d2012l) and lChen et alj 



are dis cussed in 

We treat stochastic acceleration in the whole blob as a dif- 
fusive term in the Fokker-Planck equation, while the putative first 
order Fermi acceleration at the shock front is simplified as directly 
injecting high energy particles with a power law distribution into 
the zones where the shock is present. 

The main parameters of the model are (see also Table [T}: 

• r, the bulk Lorentz factor of the jet. 

• R and Z, the radius and height of the cylindrical simulation 
region. 

• B, the magnetic field strength, assumed here to be homoge- 
neous and non variable. 

• n c , the initial electron number density. 

• t' acc and t' csc , the time-scales parameterizing the stochastic dif- 
fuse acceleration and particle escape. 

• Qpick, the rate at which the blob constantly picks up mildly 
relativistic electrons (with a narrow Gaussian distribution centered 

at 7 P ick- 

• LLj , the luminosity of the relativistic electrons injected by the 
shock. 

• Pinj, 7min,inj and 7max,inj are the spectral index, minimum 
and maximum Lorentz factor of the electrons injected locally by 
the shock with a power law spectrum. 

• ^?blr, /blr are the size of the broad line region, assumed to 
be spherical and the fraction of the luminosity from the accretion 
disk that contributes to the radiation energy density within its vol- 
ume. i?m and /ir are the corresponding parameters for the case 
of the infrared emitting torus. Collectively we also call them R cx t 
and / cxt without distinction between BLR and IR. We will discuss 
them in the following sections. 

At the beginning of a simulation each zone of the blob is filled 
with electrons with density n c and a power law spectrum with slope 
and energy range given by p, 7 m i n and 7 max . The same parameters 
are used in every zone. These parameters are not listed in Table Q] 
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because their relevance is minimal; they are simply seeding each 
zone with electrons at beginning of the simulation, and they do not 
represent the electron distribution once the simulations starts. The 
actual electrons distribution in each zone will be different from the 
seeded one, and from zone to zone, as they evolve separately. In 
the steady state it would approximate a power-law like distribution, 
with 7 m in determined by the energy of the picked-up electrons. The 
spectral index is related to t^ cc and t' csc , with faster escape and 
slower acceleration generally leading to softer spectrum. The 7 max 
is determined by the competition between acceleration and cooling, 
with faster cooling meaning smaller maximum electron energy. 

3.2 External Radiation 

The relativistic jets in Active Galactic Nuclei (AGN) may reside in 
dense external radiation environments, with contributions fr om the 
BLR (Tave cchio&Ghiselli ni 2008; P outanen & Sternl2010h or the 
warm dust of the molecular torus hypo thesized to exist beyond the 
accretion disk dMalmrose et al. 201 ll). For a thorou gh discussion 
we refer the reader tolGhisellini & Tavecchiol (l2009). 

The dominance of different sources of external radiation is 
connected to the location of the 7-ray emitting region within the 
jets. The radiation from the BLR can be dominant only when 
the emission region is located at sub-parsec distance from the 
central engine of the AGN. Beyond that distance, the infrared 
radiati on from the dust torus is likely to dominate on par sec s 
scale dGhisellini & Tavecchioll2009l) . IPoutanen & Stern! d2010l) ar- 
gue that the GeV spectral breaks of FSRQs observed by Fermi/LAT 
are a sign of 7-ray absorption inside the BLR. Meanwhile, 
iMarscher et alj d2010) used the correlation between radio knot ap- 
pearance and 7-ray flares to identify the location of the emission 
region at several parsecs from the central engine. We will test the 
viability of both of these two sources of external photons, and see 
if they can produce the SEDs and light curves observed. 

We regard the big blue bump clearly visible in the SED of 
PKS 1510—089 when in its lower brightness states as unbeamed 
thermal emission from the accretion disc, and match the data with 
a luminosity of 4 x 10 45 ergs/s and a temperature of 3 x 10 4 K. 

This disc emission is used to estimate the energy density expe- 
rienced in the blob rest frame while its location is within the radius 
of BLR, J?blr, according t o the following transformation equation 
dGhisellini & Madaulll996h : 



Um 



17 /blr Ld r 2 
12 4^i?| LR c ' 



(1) 



where La is the disc luminosity. 

Similarly, assuming for simplicity that the region where radia- 
tion from the dusty torus yields a significant energy density can be 
described as a spherical volume of radius Ljr, the energy density 
within this region as seen in the blob rest frame can be estimated 
by 

/m L d r 



ui. 



47Ti? 2 R C 



(2) 



dGhisellini & Tavecchioll2009l) . 

For the spectrum of BLR emission, we consider two cases: an 
approximation as a plain single temperature blackbody peaked at 
1.5Ti/Lya, or a more realistic spect rum obtained by taking t he un- 
beamed BLR spectrum in Fig. 4 of iTavecchio & Ghise llini (2008) 
and beam it according to the equation: 



here v\ = i/'/[r(l + /?)], 1/2 = v'/T, I(v) is the unbeamed inten- 
sity spectrum, v and v' are the frequency in the observer's frame 
and blob frame respectively. 

For the infrared emission from the hot dusty torus, we use 
a bla ckbody spectrum with temperature dGhisellini & Tavecchiol 
|2009» : 



T' IR = 370 r(l-/3cosa)K~ 370 rK, 



(4) 



u{v) =wc v L 



(3) 



where a is the angle between the jet axis and the line connecting 
the source and the jet, so a ~ n/2 for the torus. 

For computational ease, we simplify the model by assuming 
that all the external photons are traveling in the upward direction 
in the frame of the blob, as illustrated in Fig. [T] All the external 
photons enter the blob through the lower boundary and the external 
flux is just the energy density times c. This is a valid approximation 
for the typical (large) values of the Lorentz factor appropriate for 
blazars, which we adopt for ease of computation, although our code 
allows to setup the flux from the external illuminating source with 
an angular distribution. 

To produce the observed SEDs the disc emission is added as a 
non-varying component to the beamed emission from the jet in the 
post-processing of the simulation results. 



3.3 About model parameters freedom and constraints 

In the EC models, we have 5 basic observables (variability time- 
scale T va r, synchrotron luminosity L sync , estimated IC peak fre- 
quency ^ic,p, IC luminosity Lie, and 7-ray spectral index a 7 ) to 
constrain 6 free parameters (R, B, n c , / GXt , 7min and pi n j). How- 
ever, there are in fact additional constraints available to further limit 
the usable range of parameter space. For example, as we will dis- 
cuss later, in the EC cases the SSC flux can not be too high in order 
to be consistent with the moderate X-ray variability, which in turn 
translates into a requirement on the magnetic field strength to be 
sufficiently large. On the other hand, the required diffuse stochas- 
tic acceleration time-scale (t' acc ) and hence the cooling time-scale 
can not be too short, which then imposes a limit on how large the 
magnetic field can be. 

Moreover, as already noted, the value of p and 7 max in the 
quiescent state are not direct input parameters. They are the results 
of the combination of i acc and t' csc (and the relevant cooling time- 
scale). The synchrotron peak frequency f S ync,p is not always used 
as a constraint because for the physical conditions of our simula- 
tions it is not a result of the electron distribution, but it is often 
determined by synchrotron self-absorption. At the same time, ob- 
servationally its position in the SED is only poorly constrained by 
currently available data. 

For some other parameters, such as the bulk Lorentz factor of 
the jet, the radius of the BLR or torus regions, values are set on the 
basis of empirical estimates obtained by independent studies and 
therefore may not be regarded as truly free parameters. On the other 
hand, because the energy density of the external radiation in the 
blob frame, U' cy . t ~ / ex t r 2 /L! 2 xt , and L E c/L sync ~ U^ t /B 2 , a 
given SED or set of SEDs will impose a relationship between these 
parameters. 

The size-scales of BLR and torus is generally found to fol- 

1/2 

low a relationship like R ~ L d j sc , with some range in the nor- 
malization typically yielding a -Rblr = few x 10 17 cm and 
Rm — few x 10 18 cm for a disc luminosity of the order of that 
inferred in PKS 1510-089 on the basis of the observed BBB. We 
thus decided to set _R cxt to the value closed to the ones expected 
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Time in days [obs. frame] log(E) [eV] l°g(E) [eV] 

03468 10 12 14 -3 0369 -3 0369 




200 400 600 800 1000 1200 12 15 18 21 24 12 15 IB 21 24 



Time since simulation start [ks, obs. frame] log(i/) [Hz] log([/) [Hz] 

Figure 2. Left, middle: Light curves and SEDs for the quiescent state of the EC/BLR model, using a blackbody approximation for the BLR spectrum. The 
histograms in both figures are the results of our simulation, with the energy band chosen shown in the legend of the left figure. The SED snapshots times 
are shown in legend of figures and marked with matching color segments in the light curves p lot. In this and all the paper SED plots the data points are 
multiwavelength SEDs of PKS 1510-089 mostly in the spring of 2009 from lAbdo et al.1 l l2010al) and lD' Ammando et al] fed lh . The optical/infrared points, 
also show in the inset, are for the following dates: blue for March 10, orange for M arch 18, black for March 19, red for Match 25/26 (flare peak). The 7-ray, 
Fermi/LAT spectra are: blue squares for t he quiescent state of the early 2009 (see lAbdoet al. 2010a), red squares the end of March 2009 flare, black empty 
squares an earlier weaker flare (flare 'a' in AbdiLei_aL2i)L0a!!) which we show as a plausible reference for the gamma-ray state before and after the March 2009 
flare. The gray points in radio, submm and X-ray are not strictly simultaneous but we regard them as representative because of the very modest variability 
exhibited by the source during those months in these bands. The X-ray data include the XRT data avera ged during the March 20 09 flare, and five year integrated 
BAT data in hard X-rays. Right: The SEDs for the quiescent state of the BLR model using the Tavecchio & Ghisellini 12008) BLR spectrum (parameters are 
listed in TableQ] nf /blr). 



from these estimates and given the uncertainty on the / cxt to re- 
gard its value as a free parameter controlling the normalization of 
the EC emission. 

4 RESULTS 

As illustrated in Section|2]our aim is to reproduce the quiescent and 
flaring states of PKS 1510-089. 

We model a flare as being caused by an injection of relativistic 
electrons ascribed to the effect of the interaction of the blob with 
a standing shock, as described in Section |3T| As the shock travels 
through the blob (in the blob frame) it injects new particles locally, 
i.e. in the zones where it currently resides. These newly injected 
electrons are treated in the same way as all other electrons in each 
zone; they radiate through synchrotron and IC, and evolve accord- 
ing to the same Fokker-Planck equation. 

We will review cases with BLR or torus as the dominant 
sources of external photons for EC, as well as a pure SSC model, 
and compare light curves and SEDs with those of the March 2009 
observations of PKS 1510—089. In discussing the comparison we 
will mainly focus on the SEDs, in particular on the infrared, X-ray 
and 7-ray bands. The light curves, which we show for several ob- 
servable bands, provided an important constraint guiding the anal- 
ysis and identifying suitable parameter values but they do not illus- 
trate the differences between different cases as clearly as the com- 
parison of data and SEDs for multiple epochs. This is also due to 
the fact that purely in terms of intensity variation in narrow bands, 
the longest time-scale dominates (modulates) the time profile of a 
flare, and in the cases discussed here the source crossing time is 
larger than time-scales for electron processe s. This was not true for 
instance for Mrk 42 1 , for which as shown in lChen et al.l d2011bl) the 
light curve profiles for various X-ray and TeV bands were different 
depending on whether the particle related time-scales were faster or 
slower than the source crossing time, providing us with additional 
diagnostics. 



4.1 Quiescent state: no shock 

First we try to reproduce the quiescent state SED using as reference 
the period of the last three months of 2008 during which the 7-ray 
flux measured by Fenm'/LAT remained low and with small varia- 
tions ( the 'quiescent' period in the naming adopted in lAbdo et alj 
l2010j) . We begin with using the blackbody approximation for the 
spectrum of the BLR emission. We show the results of this simula- 
tion in Fig. [2] with the parameters used listed in TableQ] 

In this simulation particle injection is exclusively given by the 
steady-state pick-up term Qpick- Because there is no flaring activity, 
the flux level at every wavelength reaches the steady state after a 
few light crossing times, and the light curves remain almost flat 
except for statistical fluctuations. 

It is interesting to note that the R-band light curve reaches a 
flux level higher than the quiescent level before it reaches steady 
state (effect noticeable also in the SED in the infrared band). This 
is because the external photons need some time (of the order of 
one light crossing time) to diffuse through the whole blotQ. Dur- 
ing that time some zones which have not yet received the exter- 
nal photons will experience significantly less IC cooling (which is 
dominant over synchrotron cooling in this case). Hence the higher 
energy portion of their electron spectrum will remain at a relatively 
higher level and produce a relatively more intense synchrotron ra- 
diation during this phase. This is an example of how the light travel 
time effect can affect the actual physics in the jet, not only just 
the way we perceive the emission. While in this particular case the 
initial phase of the simulation and particle evolution is not of astro- 
physical interest, this effect is realistic and illustrates one important 

1 The same is true for the internally produced synchrotron photons scat- 
tered by SSC. Their contribution to the photon field in each location will 
need a time of the order of the light crossing time to stabilize, or, in the 
case of variable synchrotron emission, to respond to the changes of inten- 
sity happening elsewhere in the blob. 
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Figure 3. Light curves (left) and SEDs (right) for the first EC/BLR scenario. In these and all following light curve plots the blue circles in the lower panel 
show the observed R-band light curve in March 2009 from the GLAST- AGILE Support Pro gram (GASP), those in the upper panel show the Fermi/LPH light 
curve above 0.2 GeV in March 2009, simultaneous with the R-band data (Abdo et al. 2010a). The SED data points are the same as those in Fig. [2] Simulation 
light curves are dotted before T = to emphasize that that interval should be regarded as a setup time for the simulated region. The two grey vertical dashed 
lines mark the start and end of the injection time, i.e. the time during which the shock is crossing the blob. The SEDs correspond to the times (in the observer 
frame) given in the legend, and marked with colored ticks in the light curves plot. Model parameters are given in Tabled] blrl5. 



aspect of taking into account the finite size of the source and the ef- 
fect of light travel time on the physical evolution of the system and 
its emission. 

The SEDs match the low state data points at optical and 7-ray 
frequencies pretty well. The radio data points do not match because 
it is likely that the radio emission comes from additional emission 
regions rather than just the one producing the optical and 7-ray 
emission. The simulated spectrum in the X-ray regime is much 
harder than the observed one. This improves significantly by using 
a more d etailed description of the BLR spectrum such as that dis- 
cussed bv lTavecchio & Ghisellinil fe008). as shown by the SEDs in 
Fig. [2] where we used a BLR spectrum obtained from their analy- 
sis. This confirms their conclusion that an accurate treatment of the 
BLR spectrum used for the EC emission is necessary in produc- 
ing the spectrum in the soft to medium X-ray band (~0. 1-10 keV), 
where usually high quality data are available. This band is of crit- 
ical importance because it can provide constraints on the relative 
contribution of SSC and EC, and also on the characteristics of the 
lower energy end of the electron distribution. 



4.2 EC/BLR: shock crossing with T = 15 

Starting from a baseline quiescent state like the one just discussed, 
we model the flare for a case in which the dominant EC contribution 
is from the BLR. For this first case we adopt a blob bulk Lorentz 
factor r = 15. The results are shown in Fig. [3] Light curves show 
that the optical/infrared variability is larger in the R-band than in 
the bluer Swift/UVOT bands, due to the contribution of the non- 
variable big blue bump emission in the UV. This is consistent with 
observations. 

Both the X-ray light curve and the SEDs clearly show that our 



simulations predict large-amplitude variations in the X-ray flux. 
This is at odds with observations which show only modest X-ray 
variability (less than a factor of 2) throughout the entire 2009 ob- 
serving campaigns, including the largest flares. The excessive X- 
ray variation in our simulations is the result of SSC emission. Al- 
though currently our model does not track separately photons of 
EC and SSC origin, the parameters we use indicate that emission 
by SSC is not negligible in the X-ray band even in the quiescent 
state. In case of a flare caused by changes in the electron spec- 
trum and/or density the amplitude of variation of the SSC emission 
will always be larger than that of the synchrotron (as long as we 
are considering electrons scattering in the Thomson regime, which 
is the case here). Therefore as we model the factor of 10 increase 
of the non-thermal emission in the optical/infrared band the corre- 
sponding IC emission will vary by a factor up to a 100, with a large 
contribution in the X-ray band. 

This SSC variation makes this model and parameter set not 
consistent with the observed features of PKS 1510—089. This ex- 
ample also illustrates the importance of modeling the time evolu- 
tion of the SED, rather than simply modeling with sets of param- 
eters left fully free to vary between different epochs the high and 
low state SEDs, or even sequences of SEDs taken close enough in 
time during a single flare that they are very likely to be related to 
each other as part of the development of a single event. 



4.2.1 Quiescent vs. flaring state and the importance of 
time-dependent multi-zone modeling 

The comparison between the SEDs of the steady and flaring states 
shown in Figures [2] and [3] illustrates another important difference 
between a time-dependent multi-zone simulation and a one-zone 
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Figure 4. Same of Figure[3]for the second case of EC/BLR flare simulations, with T = 15 and higher minimum energy for the injected electrons, 7 m i n = 90 
(see Table[T] blrl5highgmin). 



(effectively point-like) simulation, for which the quiescent case can 
be considered a prox}0: the SEDs of the more realistic model are 
significantly more complex in the shape and variation of the high 
energy component. This is the result of both internal and external 
LTTE giving to the observer a mix of emission produced at different 
times in zones at different stages of the flare development and with 
electron distributions at different stages of their evolution. 

Even if locally the processes affecting the electrons are fast 
and the particle spectrum could be regarded as reaching rapidly a 
steady state (in case of injection lasting for a long enough time), 
the sequence of SEDs produced i n a flare is not equi valent to a 
sequence of steady state SEDs (see lBonnoli et alj|201(t for an ex- 
ample of this approach applied to the FSRQ 3C 454.3). 

Admittedly in this paper we are presenting one possible sce- 
nario for the flaring state of a FSRQ, but this type of differences 
can be expected to exist for a wide range of plausible scenarios of 
variable emission from a relativistic jet. 



4.3 EC/BLR: T = 15, with higher 7 min 

A potential remedy for the excessive X-ray variability crisis could 
be to increase the minimum energy of the injected electrons. The 
peak of the emission of the variable component, driven by the 
electron injection, is determined by the electron's 7 m m- A higher 
7min may push the SSC emission by these electrons to peak at a 
higher frequency where it can be 'hidden' beneath the rapidly ris- 
ing stronger EC emission. However, 7 m i n is fairly constrained by 
the observed MeV-GeV 7-ray spectrum, namely by the fact that we 
do not observe a spectral turn over in the FermifLXT data. Too large 
a 7min will produce a EC SED peaking in the FermifLAT observa- 



2 The steady state emission is effectively equivalent to what would be pro- 
duced by a one-zone non-time-dependent code. 



tional band. With this in mind, we tested one case with a slightly 
higher injected 7 m i n = 90. The results are shown in Fig. [4] 

The peak of the EC SED in this case is shifted close to the 
lowest energy data point of the FermifLAT spectrum, indicating 
that 7 m in = 90 is already of the order of the largest value that 
we can use with the current BLR spectrum and a Lorentz factor of 
15. On the other hand, while there is some change on the spectral 
shape of the high state X-ray spectra, the fundamental problem of 
the excessive amplitude of its variation is not mitigated. The in- 
consistency between the observed and simulated X-ray variability 
remains a problem in the higher injected 7 m i n case. 

4.4 EC/BLR: higher Lorentz Factor, T = 25 

Instead of moving the SSC emission in frequency, another route 
to solve the X-ray variability inconsistency may be to decrease the 
level of the SSC emission relative to the other components. In or- 
der to do this while keeping the same level of synchrotron emission, 
we need to decrease the ratio between the synchrotron (SSC seed 
photons) energy density and the magnetic energy density. We can 
achieve this by increasing the Doppler factor because in that case 
the synchrotron energy density in the blob frame needs to decrease 
accordingly to match the optical data. This requires to change the 
values of other parameters to produce a SED well fitting the ob- 
served one. The modified parameters are reported in Table [TJ The 
results are shown in Fig. [5] 

In this case, the luminosity in the X-ray dip between the two 
main components is indeed lower. Nevertheless because any SSC 
emission occurring in this band would still be varying by a factor 
larger than that of the optical flux, the range of the X-ray variation 
remains large, and easily exceeding the constraint set by the well 
measured X-ray intensity and spectrum. However, since in this case 
the quiescent X-ray flux is lower than the observed one, it may 
be possible to explain the X-ray band spectrum as comprising a 
contribution from additional, relatively cooled blobs, which would 
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Figure 5. Light curves (left) and SEDs (right) for the EC/BLR case with higher bulk Lorentz factor, F = 25. Parameters are given in Table[TJ blr25. 



partially dilute the large variation. However, even taking that into 
account, considering the spectral variability present in the varying 
X-ray component, it may not be straightforward to reconcile the 
overall X-ray properties with the remarkably stable observed spec- 
tra. 

4.5 EC dominated by IR emission from the torus 

Emission from hot dusty, molecular, gas in the putative torus sur- 
rounding the accretion disk is another plausible source of exter- 
nal photons in the immediate environment of the relativistic jet. 
This scenario is motivated by the observation of the coincidence 
of 7-ray flares wit h the appearance of ne w knots in radio images 
of PKS 1510-089 tlMarscher et"aT]|20irj|) . These observations hint 
that the emission region responsible for 7-ray flares is located at 
parsec scales, which is beyond the usuall y inferred radius of the 
BLR dKaspi et alJl2007l;lBentz 
the IR torus jPier & Kroliklll992bl 



et alf: 
J2b.a) 



usually mien 
200rj,l2009). 



At this distance, 
becomes the main candidate 
as the source of E C seed photons (Ghisellini & Tavecchioll2009l : 
ISikora et al.ll2009l) . 

We calculate the energy density and temperature of the torus 
emission according to equations ((2} and Because the energy 
density and SEDs of the emission from the torus are very different 
from those of the BLR, it is necessary to change the values of sev- 
eral parameters to be able to match the quiescent and then flaring 
states. The best set of parameters is reported in Table [T] and the 
light curves and SEDs in Fig. [6] 

The broadband SEDs compare reasonably well with observa- 
tions; the light curves vary on a time-scale consistent with the data; 
the optical and 7-ray light curves are well correlated with no signifi- 
cant lags; the variations in the optical and 7-ray bands have similar 
amplitude; the variations in the UV band are less prominent than 
those in the optical. However, it is worth noting that it has proven 
to be difficult to concurrently match the GeV spectrum and the soft 
slope through the infrared bands. Although we can not claim to 
have achieved a perfect match in the latter in the previous cases, 
this is an issue that did not seem to emerge as seriously as here. 



The problem of X-ray variability in this case is similar to the 
one in the BLR case. It is likely that similarly to the previous case, 
a higher Doppler factor would lower the flux produced by SSC in 
the X-ray band. 

One of the main differences between the torus emission and 
the BLR emission as source of EC seed photons is that the ow- 
ing to its lower temperature the torus radiates at lower frequency 
compared to the BLR. This means that to scatter these seed pho- 
tons to the same 7-ray energies, the energy of the electrons needs 
to be higher than those needed in a BLR scenario. In the context of 
the scenario presented in this work, this means that more efficient 
particle acceleration is needed to sustain the high energy electrons, 
which have faster cooling times. This also means that faster parti- 
cle escape is to be expected in the torus scenario, otherwise the ac- 
celerated particles will not form a power-law distribution that can 
produce emission with the observed spectral shape. It turns out that 
the value for the particle escape time-scale parameter needed in this 
case may be too fast to be realistic (t' CBC = 0.015 Z/c). However, it 
is worth emphasizing that, chosen name notwithstanding, the 'es- 
cape' term in the kinetic equation may be regarded as a crude ap- 
proximation for a describing a generic energy independent process 
affecting electrons, and in this sense a value significantly smaller 
than the source crossing time may not be automatically be consid- 
ered unphysical. Nevertheless, because the time-scale for plausi- 
ble candidate such as actual escape or adiabatic cooling processes 
would likely be of the order of Z/c, it would certainly be desirable 
to not be forced to such extreme values for t' eBC , and we regard this 
as a serious problem for EC on torus emission, in the framework 
discussed here. 

4.6 Pure SSC 

We already cited some of the recent results suggesting that 
the active, 7-ray emitting region, in the jets of several blazars 
may be located beyond the size-scale of the BLR on the ba- 
sis of correlated multiwavelength variability and VLBA imag- 
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Figure 6. Light curves (left) and SEDs (right) for the case with the EC dominated by radiation from the torus. Parameters are listed in Table[T] torus 1 5. 
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challenges traditional EC scenarios because jet emitted TeV pho- 
tons would be readily lost by photon-photon pair production with 
the copious soft photons surrounding the jet. These findings have 
thus stimulated a renewed interest in the possibility that even for 
some of these powerful jets in systems with luminous accretion disk 
and broad line emission, the 7-ray emissi on may be predomina ntly 
by SSC (for analysis on 3C 279 see e.g.. lB6ttcher et al.ll2009l) . Fi- 
nally, studying a large sample of well characterized blazars detected 
by Fermi/LAT and with an estimate of their jet intrinsic power, 
lMeveretalJj2012l) find that, while for the highest jet power objects 
there is a clear collective sign of EC being the dominant mechanism 
for their 7-ray emission, the properties of the rest of the population 
are consistent with SSC. 

Moreover, as discussed in the previous sections, our simula- 
tions of EC models for the type of flaring blob scenario presented 
here show large variability in X-ray which is inconsistent with the 
very robust observational finding of lack of significant variations 
even during the flaring phase, and which is caused by the SSC con- 
tribution. The existence of this latter is unavoidable, and our tests 
suggest that even mitigating its effect is not a trivial endeavor (§ 14.31 

§03). 

Therefore, it seems natural to want to test if we could repro- 
duce the benchmark observations with a pure SSC scenario. The 
results are shown in Fig. [7] 

The SSC model reproduces well several aspects of the obser- 
vations. In particular, the amplitude of the flux variability and the 
spectral change in the X-ray band are overall much closer to what 
was observed. The Fermi/LAT band spectrum and intensity is also 
well matched. Because of the presence of the steady disk emis- 
sion mostly contributing to the blue/UV flux, the flux in the R-band 
varies less than the SwiftfUVCfT band. 



However, the simulations produce an X-ray spectrum consis- 
tently harder than the one observed, mostly determined by the sharp 
'cut-off of the lower energy end of the electron distribution. The 
observed SED, namely the very hard X-ray spectrum, constrains 
this latter to a fairly high 7 m in, much larger than the values used 
in the EC cases. The simulated infrared spectrum is also somewhat 
harder than the spectrum observed in the intermediate state. Finally, 
the frequency at which the synchrotron spectrum peaks tends to be 
too high, which is a consequence of the high 7 m i n required by the 
X-ray spectrum. Since the SSC model has fewer free parameters 
than the EC model, they are more constrained than those in the EC 
models, and do not leave us much freedom for improving signifi- 
cantly on these issues. 



5 DISCUSSION AND CONCLUSIONS 

We have modeled multiwavelength variability produced in the jet of 
PKS 1510—089, with EC model involving external radiation from 
BLR or dusty molecular to rus, as well as pure SSC model . 

As discussed also by iTavecchio & GhisellinJ d2008t) . the use 
of a realistic BLR spectrum seems indeed to be critical to model 
accurately the inverse Compton component, in particular its lower 
energy end which is observed in great detail in the X-ray band, 
thus providing a powerful diagnostic on model parameters. The 
results presented here, namely the challenging issue of the large 
X-ray variability caused by the SSC contribution, which has a dis- 
tinctly different spectral shape from the EC, further strengthen the 
importance of modeling as accurately as possible the source of the 
EC seed photons. 

In the same context, we should comment on the fact, noted 
in Section [331 that in order to obtain reasonable SEDs and light 
curves in the EC/BLR cases we had to adopt a reprocessed frac- 
tion of BLR luminosity, /blr, that is significantly smaller than 
what one would expect . A more traditional /blr, such as 0.1 in 
iGhisellini & Tavecchid J20091) . would yield a more intense BLR 
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Figure 7. Light curves (left) and SEDs (right) for the SSC case. Parameters are listed in Table[TJ ssc. 



radiation energy density and the increased electron energy loss rate 
in the blob would push the model towards fairly extreme values 
for several of the parameters. This could be regarded as an indica- 
tion that the general scenario that we studied in this work may not 
be a suitable explanation for flares in FSRQs. On the other hand, 
this scenario yields acceptable, though not perfect, results for the 
EC/torus and the SSC cases. 

As discussed in Section 13.31 the combination of 
(/ext r 2 /7?ext)/-B 2 is related to the ratio between the IC 
and the synchrotron luminosities (Compton dominance), which is 
an observable. Although its value may vary during a flare, the value 
of the Compton dominance imposes a relationship between those 
parameters of the source (hence model). A given source or source 
state may have a preferred value for the Compton dominance 
which should be matched by a combination of / cxt , V and i? oxt 
(and B), independently on whether one wants to model it with 
EC/BLR or EC/torus. For the Compton dominance ratio that seems 
to work well for PKS 1510—089 it is possible to adopt plausible 
values or i? cxt and / cxt for the EC/torus case and less so for the 
EC/BLR case. This consideration does not depend strongly on 
the scenario adopted for the the blob and the mechanism causing 
outbursts and it may suggest that indeed PKS 1510—089 may not 
be modeled by EC/BLR in general. Alternatively, it could mean 
that our naive idea of the BLR and torus does not represent their 
real structure and causes the EC/BLR model to fail to match the 
data. Another possibility is that, the relativistic jet in this source 
has a milder bulk Lorentz factor F than 15 or 25 that we adopted 
in this paper. However this in conflict with the high T required to 
avoid the large SSC contribution. 

Leaving aside the above considerations, and looking at the re- 
sults of the simulations as run with the best fitting parameter sets, it 
is difficult to identify a clear superiority for any of them. None can 
be convincingly excluded. 

The EC/BLR scenario suffers the problem of the large X-ray 
variability, and it produces better, marginally satisfactory, results if 
we adopt a larger value for the bulk Lorentz factor, which is still 



within the observed range. It also requires parameters for the BLR 
which are at least unusual. 

The EC/torus also has the problem of highly variable SSC in 
the X-ray band, which we can expect to be mitigated by using a 
larger Lorentz factor, and it works with reasonable parameters for 
the external radiation component, but it requires very fast parti- 
cle acceleration and escape time-scales to maintain the balance be- 
tween acceleration and cooling for the high energy electrons neces- 
sary to scatter to the GeV band the lower energy external photons. 

The pure SSC model naturally addresses the X-ray variability 
issue, but at price of very high minimum electron energy to avoid 
significant emission in X-rays, which in turn make it difficult to 
reproduce the synchrotron peak region SED and makes the X-ray 
spectrum significantly harder than the observed one. Because its 
GeV emission comes from scattering lower energy photons than 
the EC/BLR case, like the EC/torus case the SSC requires fast es- 
cape time-scale, though less extreme, which could be considered a 
serious problem. 

In all three cases, their X-ray problems may be relaxed if there 
is a slowly varying contribution from other emission regions, filling 
the dip between the synchrotron and IC components, softening the 
X-ray spectrum and diluting its variations. 

Overall, while within the framework studied in this work we 
can not firmly discriminate the three main scenarios, and the frame- 
work that we adopted is just one possibility, our results show clearly 
the differences produced by a more realistic treatment of the emit- 
ting source in the shape of SEDs and their time variability over 
relevant, observable time-scales. 

These results demonstrate that proper modeling of the high 
quality data produced by the plethora of best multiwavelength cam- 
paigns on the brightest FermifLKT blazars to exploit the wealth of 
information that they carry and advance our understanding of their 
physics can only be achieved with time-dependent multi-zone sim- 
ulations. 

Looking forward, there are several aspects of the predicted 
source variability that we have not discussed, such as more de- 
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tailed analysis of correlated variability (e.g. time lags). Moreover 
as shown here, SEDs exhibit complex features and variations and 
it may be possible to explain the interesting, challenging, obser- 
vation of large break in the FermifLKI spectra seen in a hand- 
ful of sources, com prising both FSRQs and low-peaked BLLacs 
dAbdo et all boiOch . The observed breaks are larger than what 
would be produced by cooling and too sharp to be consistent 
with an exponential cutoff. Several ideas have been proposed to 
explaining them a s due to external factors, namely 77 absorp- 
tion outside the je t jPoutanen & Sternl20ld : lFinke & Dermerl20ld: 
lAckermann et al.l|2010h . and some are somewhat inconsistent with 
other multiwavelength properties. 

Finally, in these simulations we have used a simple model of 
energy independent particle acceleration, whereby in order to main- 
tain a power-law electron distribution with a certain slope, the parti- 
cle escape time-scale is tied to the acceleration time-scale. The fast 
radiative cooling in our model requires fast particle acceleration 
and hence fast particle escape. This required escape time-scale has 
turned out to be smaller than the limit of Z/ c, which happens if all 
particles travel freely at the speed of light. However, the actual par- 
ticle acceleration process is expected to have a more complicated 
energy dependence than assumed here. Therefore, the acceleration 
and escape time-scales derived here are only for instructive pur- 
poses. We will investigate whether a more realistic energy depen- 
dence of the particle acceleration process may change the required 
acceleration and escape time-scales. 
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